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Abstract 

A  delta  wing  at  a  high  angle  of  attack  produces  two  vortices  that  generally  undergo 
dramatic  increases  in  core  size,  followed  by  the  formation  of  regions  of  reversed  flow  .  Th:s 
phenomenon  is  called  vortex  breakdown  and  can  have  significant  effects  on  the  aircraft's 
lift.  drag,  and  moment  coefficients.  The  objective  of  this  thesis  is  to  provide  a  baseline 
model  of  the  compressible  vortex,  independent  of  the  complex  body  interaction  with  the 

:  r.  -  i.  •.  * 

delta  wing.  The^model  is  then  used  to  simulate  vortex  breakdown  for  various  vortex 
strengths,  Reynolds  numbers,  and  Mach  numbers  with  particular  attention  given  to  the 
effects  of  compressibility. 

After  running  many  simulations  it  was  found  that  Mach  number  has  a  favorable  effect 
by  delaying  vortex  breakdown  as  defined  above.  Holding  Reynolds  number  and  vortex 
strength  constant  while  increasing  Mach  number  reduced  the  effective  vortex  strength 
while  compressing  the  flow.  Another  important  result  of  this  compressible  flow  study  was 
the  disappearance  of  non-unique  solutions  at  Re  =  200  and  V  =  1.0  as  Mach  number  was 
increased.  No  paths  of  non-unique  solutions  were  found  for  A/  >  0.2. 

-V  ■  -  -  ,  ^  ",  r>  .  *  ■  f  > 

,  *’  ‘  / 
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NUMERICAL  SIMULATION  OF  COMPRESSIBLE  VORTICES 

/.  Introduction 

Vortex  breakdown  is  a  phenomenon  that  occurs  in  rotational  flows  with  a  concen¬ 
trated  core  of  vorticity  within  a  primarily  irrotational  flow.  Swirling  flows  in  pipes  and 
leading-edge  vortices  over  delta  wings  are  examples  of  this  flow.  The  term  “vortex  break¬ 
down"  is  defined  by  Beran  (2)  as  rotational  flow  with  the  development  of  a  stagnation 
point  on  the  core  of  a  vortex,  followed  by  a  region  of  reversed  flow  and  a  dramatic  increase 
in  core  size. 

Vortex  breakdown  of  the  leading  edge  vortices  over  a  delta  wing  can  have  significant 
effects  on  the  aircraft’s  lift,  drag,  and  moment  coefficients.  For  this  reason  breakdown  has 
been  studied  by  many  aerodynamicists  since  it  was  first  discovered  in  1957  by  Peckham 
and  Atkinson  (33).  A  brief  historical  background  and  an  outline  of  the  study  is  contained 
in  this  chapter. 

1.1  Historical  Background 

In  the  past  thirty-two  years  since  the  discovery  of  vortex  breakdown,  there  have 
been  many  experimental,  analytical,  and  numerical  investigations  of  the  vortex  breakdow  n 
problem.  Several  review  papers  have  been  written  on  the  subject,  most  notably  the  reports 
by  Hall  (17),  Leibovich  (26),  Leibovich  (28).  and  Hall  (14).  After  Peckham  and  Atkinson's 
discovery.  Elle  (5)  and  Lambourne  and  Bryer  (25)  conducted  further  experimental  inves¬ 
tigations  of  leading-edge  vortices.  Useful  information  was  presented  in  these  studies  but 
quantitative  results  were  difficult  to  obtain  because  of  the  sensitivity  of  the  vortex  core  to 
probe  disturbances. 

Qualitative  results  were  obtained  by  visualization,  which  led  them  to  the  discovery 
of  two  distinct  types  of  vortex  breakdown.  In  Hall’s  review  (17),  Figure  2  reproduces  the 
picture  taken  by  Lambourne  and  Bryer  of  vortex  breakdown  in  the  leading-edge  vortices 
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over  a  delta  wing.  Flow  visuaJization  was  made  possible  by  injecting  dye  into  the  fluid  at 
the  wing  apex.  With  this  technique  they  discovered  “bubble  vortex  breakdown"  which  is 
characterized  by  the  near  axisymmetric  swelling  of  the  core  into  the  shape  of  a  bubble, 
followed  by  turbulent  disintegration  of  the  vortex.  The  other  type  of  breakdown  found 
was  “helical  vortex  breakdown"  which  can  be  described  as  an  abrupt  transformation  of 
the  linear  core  into  a  helical  filament.  Helical  vortex  breakdown  occurs  after  several  turns 
of  the  helical  core. 

The  discovery  of  bubble  vortex  breakdown  led  several  investigators  to  conduct  swirling 
pipe  flow  experiments.  The  pipes  were  made  of  Plexiglass  to  permit  flow  visualization. 
The  vortex  strength  and  Reynolds  number  associated  with  the  flow  could  be  controlled 
in  this  apparatus  by  moving  swirl  vanes  at  the  inflow  and  changing  the  mass  flow  rate 
respectively.  This  type  of  experiment  was  conducted  by  Harvey  (18),  Kirkpatrick  (22), 
Sarpkava  (35),  Sarpkaya  (36),  Sarpkaya  (37),  Faler  and  Leibovich  (6),  Garg  (8),  and  FaJer 
and  Leibovich  (7).  Quantitative  data  was  obtained  by  the  last  three  through  the  use  of 
laser  doppler  anemometry.  By  constraining  the  flow  entering  the  test  section  to  rotational 
symmetry,  Faler  and  Leibovich  (6)  were  able  to  classify  seven  distinct  types  of  breakdown 
with  five  of  them  containing  a  stagnation  point  on  the  core,  satisfying  the  earlier  definition 
of  vortex  breakdown.  The  bubble  type  of  breakdown  was  classified  type  0,  whereas  the 
helical  form  was  classified  a  type  2  breakdown. 

The  bubble  form  of  breakdown  was  the  most  attractive  type  to  simulate  numerically 
due  to  its  axisymmetric  nature.  It  could  be  modeled  in  2-D  with  a  square  computational 
grid,  saving  on  computer  resources.  Kopecky  and  Torence  (23)  were  the  first  to  model 
vortex  breakdown  with  the  Navier-Stokes  equations.  Since  1973.  results  have  been  reported 
by  Grabowski  (11),  Krause  et  al.  (24),  Nakamura  et  al.  (31),  Nakamura  et  al.  (32),  Hafez 
et  al.  ( 12),  Beran  (1),  Hafez  et  al.  (13),  Brown  and  Lopez  (3),  Lopez  (29).  and  Menne  (30). 
Nakamura  used  the  vortex  filament  method  which  allowed  him  to  generalize  the  flow  from 
rotational  symmetry.  The  other  papers  assumed  rotational  symmetry  in  the  model  of 
the  vortex  breakdown.  The  papers  by  Brown  and  Lopez  (3)  and  Lopez  (29)  compared 
computational  results  with  experimental  flow  visualizations  and  found  excellent  agreement. 
All  of  these  simulations  have  been  accomplished  using  the  assumption  of  incompressibility. 


A  compressible  study  of  the  vortex  breakdown  is  a  logicaJ  next  step  and  is  the  subject  of 
the  present  work. 

Problem  Statement 

There  are  three  main  objectives  of  this  study.  The  first  is  to  extend  the  analysis 
of  vortex  breakdown  into  the  regime  of  compressible  flow,  and  to  determine  the  trends 
of  increasing  Mach  number.  Since  delta  wing  aircraft  will  certainly  be  operating  in  the 
compressible  flow  regime,  a  compressible  study  of  the  vortex  breakdown,  independent  of 
the  complex  delta  wing  body  interaction,  is  needed.  This  objective  will  be  accomplished 
by  obtaining  a  solution  at  a  low  freestream  Mach  number.  M  ,  with  a  particular  set  of 
Reynolds  number.  Re  ,  and  vortex  strength,  V  ,  for  which  incompressible  flow  data  is 
available  and  the  vortex  breakdown  phenomenon  occurs.  A  comparison  is  made  with  this 
data  to  validate  the  code,  then  calculations  are  made  to  obtain  new  solutions  at  increasing 
Mach  number  for  the  same  values  of  Re  and  V. 

The  second  objective  is  to  investigate  the  effects  of  compressibility  on  the  occurrence 
of  non-unique  solutions  found  by  Beran  (1).  This  objective  will  be  accomplished  by  using 
Euler  Newton  pseudo- arclength  continuation  to  map  out  the  solution  path  as  a  function 
of  vortex  strength.  This  algorithm  was  developed  by  Keller  (20)  and  is  described  in  detail 
in  Chapter  III. 

The  third  objective  is  to  perform  a  grid  sensitivity  study.  This  objective  is  very- 
important  to  investigators  doing  3-D  delta  wing  simulations.  It  will  give  these  investigators 
information  on  the  grid  necessary  to  capture  the  vortex  breakdown  phenomenon  accurately. 

1.3  Study  Contributions 

This  was  a  successful  investigation  because  it  outlined  the  effects  of  compressibility- 
on  axisymmetric  vortex  breakdown,  not  previously  found  in  a  body  independent  study. 
After  running  many  simulations,  it  was  found  that  Mach  number  has  a  favorable  effect 
by  destroying  the  vortex  breakdown  phenomenon.  Holding  Reynolds  number  and  vor¬ 
tex  strength  constant,  while  increasing  Mach  number,  eliminated  vortex  breakdown,  while 
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compressing  the  fluid.  Another  important  result  of  this  compressible  flow  study  was  the 
disappearance  of  non  unique  solutions  as  Mach  number  was  increased.  Non-unique  solu¬ 
tions  were  reported  in  a  previous  incompressible  flow  study  by  Bcran  (1)  at  Re  —  200  and 
\  =  1.0.  No  paths  of  non-unique  solutions  were  found  in  the  compressible  flow  regime 
( M  >  0.3). 
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II.  Mathematical  Model 


In  this  study  of  the  vortex  breakdown,  it  was  necessary  to  use  the  full  Navier-Stoke.- 
(NS)  model  to  capture  the  compressible,  viscous  nature  of  the  trailing  vortex.  The  model 
assumes  subsonic,  compressible  flow,  and  thus  it  is  not  necessary  to  keep  the  convective 
terms  in  conservative  form.  The  NS  model  consists  of  five  equations  and  five  unknowns: 
continuity,  radial  momentum,  swirl  momentum,  axial  momentum,  and  the  energy  equation. 
The  equation  of  state  was  used  to  convert  all  other  variables  into  the  five  unknowns: 
density,  p  ,  radial  velocity,  u  ,  swirl  velocity,  r  ,  axial  velocity,  w  ,  and  internal  energy,  t  . 


2.1  Model  Assumptions 

There  were  various  assumptions  used  in  the  derivation  of  the  NS  model.  Steady-state 
solutions  were  the  object  of  this  study  and  so  ^  terms  were  set  to  zero.  To  simplify  the 
equation  of  state,  the  perfect  gas  assumption  was  used  with  constant  specific  heat  ratio 
y  =  1.4.  This  assumption  is  not  a  restrictive  one  because  vor*“x  breakdown  has  been 
observed  experimentally  at  subsonic  speeds.  This  argument  a)  “sumption 

of  Stoke’s  hypothesis,  which  is  also  good  for  subsonic  speeds. 


The  temperature- viscosity  relationship  was  derived  by  a  linearization  of  Sutherland’s 
formula: 


H(T)  = 


CiT3/3 

c2  +  r 


(2.1) 


where  for  air  at  moderate  temperatures  C2  =  110.4°/T  and  C]  need  not  be  specified  due  to 
non-dimensionalization.  First  Sutherland's  law  was  non-dimensionalized  by  p(7o).  where 
To  is  the  temperature  in  the  freestream: 


•  =  /fill  _  CiT3/?  (C2  +  Tp)  =  /  T\3/J  (Cj/Tp+l) 

*  n(To)  (C2  +  T)  cX/2  ( cyTo  +  r/To )• 

Now .  let  T*  =  T /To  and  Cj  =  C2/Tq,  to  obtain  the  non-dimensional  Sutherland's  formula. 


(2.3) 
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The  solution  algorithm  requires  that  ail  variables  appear  in  the  numerator  of  the 
algebraic  equations  resulting  from  discretization  of  the  NS  model.  For  this  reason,  a 
linearization  of  equation  2.3  is  obtained  to  ensure  that  T  appears  only  in  the  numerator. 
A  first  order  Taylor  series  expansion  was  used  in  the  linearization  process: 


•  >  /  T**  \  <  d/i 

M  =  P  (T0)  +  — 


AT*  +  O(AT’)2. 


where 


=  [£l  it  11  +  T‘ )  -  (T')3/2 

dT •  (C2*  +  T*)2  [2(y  ’  1  2  +  1  ’  ’ 


(2.4) 


( 2.5) 


but  To  =  ft  =  1  and  AT’  =  T‘  -  Tf  =  T‘  -  1. 


so 


dfT 

dT- 


3  1  1 

.2  '  (CJ  +  1)J 


(2.6) 


which  gives  the  non-dimensional  linearized  viscosity-temperature  relationship: 


2  (CJ  +  1) 


(r  - 1)+  i. 


(2.7) 


A  plot  of  equations  2.3  and  2.7  is  given  in  Figure  2.1  to  show  the  validity  of  this  assumption. 
An  analysis  of  the  maximum  difference  in  the  linear  pm  relationship  and  Sutherland’s  p’ 
relationship  was  accomplished.  At  M  =  0.4  the  maximum  error  was  0.33%  and  at  M  =  0.8 
a  maximum  error  of  2.0%  was  found. 

The  non-dimensional  thermal  conductivity-temperature  relationship  is  derived  from 
the  Prandtl  number  relationship: 

c„u 

(2.8) 


Pt=CJ£ 

k 


where  cp  and  Pr  are  assumed  constant  for  this  study.  Non-dimensionalization  of  p  and 
k  by  fio  and  k0  gives 

c„umu n 

(2.9) 


Pr  -  Cp/i 

kmk0  ’ 


but  also  equals  Pr  so 


Pr  =  %P  r. 


(2.10) 
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T* 


Figure  2.1.  Non-dimensional  Sutherland’s  Formula.  To  =  298° A' 
giving  k'  -  fim. 

Axisymmetry  is  also  an  assumption  used  in  the  governing  equations.  In  experimental 
work  two  kinds  of  vortex  breakdown  have  been  observed:  bubble  vortex  breakdown,  and 
helical  vortex  breakdown.  The  bubble  form  is  nearly  axisymmetric.  whereas  the  helical 
form  is  a  fully  3-D  phenomenon.  In  the  interest  of  geometric  simplicity,  the  bubble  form 
is  the  subject  of  this  thesis.  The  governing  equations  are  cast  in  cylindrical  form  and  then 
all  partial  derivatives  with  respect  to  0  are  set  to  zero. 

The  last  main  assumption  is  laminar  flow.  No  attempt  was  made  at  modeling  turbu¬ 
lence.  since  this  would  add  complications  to  an  already  complex  flow  field.  Vortex  break¬ 
down  has  been  observed  in  laminar*  flow  by  previous  investigators.  In  all  numerical  simu¬ 
lations  Reynolds  number  was  kept  at  200,  which  is  assumed  to  be  in  the  laminar  regime. 
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2.2  Governing  Equations 

In  this  section  the  non-dimensional,  governing  equations  are  presented.  The  deriva¬ 
tion  of  these  equations  is  accomplished  using  the  assumptions  from  the  previous  section 
and  can  be  found  in  Appendix  A.  The  non-dimensionalized  lengths  were  scaled  bv  the 
radius  of  the  vortex  core  at  the  inflow  boundary.  So  ,  and  velocities  were  scaled  by  the 
freestream  velocity.  14’  .  Variables  with  subscripts  r  or:  denote  partial  derivatives  with 
respect  to  the  independent  variables.  Following  is  a  list  of  the  five  equations  solved  nu¬ 
merically. 

Continuity: 

pu 

(pu)r  +  —  +  (pu),  =  0  (2.11  ' 

r 

u  Momentum: 

pi’2 

PUUT - —  +  p\VUt  +  (l  -  1  )(pe)r  = 

—  [cier(4/3ur  -  2/3  ^  -  2/3t it,)  +  +  u-, )  (2.12) 

+(cje  +  c2 )( 4/3tarr  +  4/3-^-  -  4/3  ^  -f  utt  +  l/3u.-,. )] 


v  Momentum: 


tr  Momentum: 


puv 

puvr  +  — —  +  pwv,  - 

1  .  V 

—  [C]er(vr  -  -)  +  CieI(t’l) 

•f  (CjC  +  CjKiVt  +  - - X  +  t’j,  )] 

r  r * 


(2.131 


pUU'r  +  pU’U'j  +  h  ~  1)(P€  ),  = 

1  ,  U 

—  [cier(uv  +  11,)  +  C(f,(- 2/3ur  -  2/3-  +  4/3u-J  (2 .14) 

+(cjf  +■  c2 )( l/3urt  +  1/3—  -f  u'TT  +  —  +  4/3u’„ )] 

r  r 
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Energy: 


'i~  +  )r  +  (Pu’e)*]  +  y  [u3  +  “i'2  +  uu'2]  +  y  [u2u'  +  r2u‘  +  u-3) 


-r^[(  u3)r  +  (ur2)r  +  (uu-2)r  +  (u2ir),  +  (v  2w).  +  (w3),  +  —  +  —  + 


U3  U  V2  UU’2, 


1  u2  V 2 

•7pkl<?r(^/3tiZ<r  -  2/3 - 2/3 UWZ  +  vvr - ■+  UtU'  -+  IT  IT,  ) 

nf  r  r 

+  Cie^uu*  +  uwr  +  t’t,  -  2/3 urw  -  2/3—  +  4/3u.’u>,) 

+  (cie  +  c 2)(4/3u?  +  4/3uur,  -  4/3urtux  +  l/3utur,  +  l/3urzu,-  -f  u2 

VVT  7  ,  7 

+  t'tVr  -  —  +  2u2tiV  +  U’r  +  U'^’rr  +  U,  +  UU„  +  V,  +  t>Vlr 

+  1/3^  -4/3^  +  4/3u-?  +  4/3u'u>„-2/3-  -  -  +  — ) 


!cje  +  c2)(err  +  Cjj  +  y  H"  cte?  +  cf 


(2.15) 


where  i?c  -  ^a,  M  -  /*»"  =  cp^,  *nd  Ci  and  c2  are  the  constants  obtained  i 

linearization  of  Sutherland’s  formula  and  conversion  of  T*  to  e: 


in 


ci  = 


r  ] 
3/2-- 


c;  +  i 


7(7~  1)M5 


(2.16) 


and 


Cl  =  1.[3/2  >  ■. 


(2.17) 


Equations  2.11-  2.15  are  then  converted  to  a  set  of  non-linear  algebraic  equations  by 
means  of  second-order  accurate  central  difference  operators,  and  utilizing  a  constant  grid 
spacing  in  both  directions: 

uI+1  -  U,_j 


fxu  = 


6  u  = 


2  Ax 

ul+]  -  2u,  +  u,_. 


(Ax)2 

On  the  boundaries,  second-order  one-sided  derivatives  were  used: 


(2. IS) 
(2.19) 


Szu  = 


3u,  -  4u,_i  +  u,_2 
2  Ax 


(2.20) 


Since  an  approximation  of  the  location  of  the  outflow  boundary  is  used,  a  second-order 
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Figure  2.2.  Swirl  velocity  typical  of  trailing  vortices. Beran  (2) 
accurate  difference  at  /  -  i  can  be  specified: 


6zu  = 


u, 


-  u.-i 

Ax 


(2.21) 


An  example  of  the  discretization  process  is  given  in  Appendix  B  by  applying  these  operators 
to  the  continuity  equation.  Once  the  governing  equations  are  discretized  the  boundary 
conditions  need  to  be  specified  to  complete  the  set. 


2.3  Boundary  Conditions 

The  flow  was  assumed  to  be  generated  by  a  delta  wing  at  a  high  angle  of  attack.  The 
resulting  trailing  vortex  is  assumed  to  produce  a  swirl  velocity  profile  downstream  of  the 
body  approximated  by  Figure  2.2.  This  profile  can  be  stated  mathematically  as 


v  = 


l>(2  -  r2)  r  <  1 
V/r  r  >  1 


(2.22) 


V  becomes  a  free  parameter  of  the  model  and  is  referred  to  as  the  vortex  strength.  When 
1=1.  the  swirl  velocity  at  the  edge  of  the  vortex  core  is  equal  to  the  freestream  velocity 
U  .  The  inflow  boundary  uses  this  swirl  velocity  profile,  uniform  axial  velocity  equal  to 
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Figure  2.3.  Characterization  of  the  flow,  Beran  (2) 

It',  and  radial  velocity  assumed  to  be  zero.  The  physical  domain  is  depicted  in  Figure  2.3 
and  the  computational  domain  in  Figure  2.4. 

The  radial  boundary  is  assumed  to  be  far  enough  away  from  the  vortex  axis  for  the 
axial  velocity  to  return  to  freestream  conditions.  Also,  swirl  velocity  is  specified  to  be 
»  V/R  and  the  radial  velocity  is  specified  in  such  a  way  as  to  satisfy  continuity  on  the  radial 
boundary.  The  outflow  boundary  is  assumed  to  be  far  enough  downstream  for  ^  =  0  for 
all  of  the  dependent  variables.  The  last  boundary  is  the  axis  of  symmetry  and  has  velocity 
components  u  and  v  equal  to  zero,  with  a  symmetry  condition  imposed  on  u\  (%£  =  0). 

The  conditions  on  p  and  e  have  not  been  described  up  to  this  point  because  different 
conditions  were  used  in  two  phases  of  the  investigation.  In  the  first  phase,  the  energy 
equation  is  replaced  with  p  -  1  to  simulate  an  incompressible  flow.  This  causes  c  to 
become  a  pressure  term  by  the  non-dimensional  equation  of  state  p  =  (7  -  l)pe.  So  the 
boundary  conditions  on  e  become  pressure  boundary  conditions.  On  the  inflow  boundary. 
Si.  (  is  determined  by  satisfying  the  u  momentum  equation  on  the  boundary.  On  the  radial 
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Figure  2.4.  Geometry  of  the  computational  domain.  Beran  (2) 

boundary  S2,  e  =  0  which  gives  p=0  (freestream  pressure).  Outflow  boundary  S3  has  the 
condition  §7  =  0.  The  symmetry  condition  §*  =  0  is  imposed  on  the  axial  boundary  S4. 

In  the  second  phase,  compressible  flow,  without  the  constant  density  assumption, 
was  modeled.  The  compressible,  columnar-flow  equations  were  solved  to  determine  an 
approximate  density  profile.  Internal  energy  is  not  specified  explicitly,  instead  the  condition 
=  0  is  used  to  enforce  a  peak  at  the  inflow.  On  the  S2  boundary,  p  =  1  is  again  specified, 
and  a  Dirichlet  condition  on  e  is  derived  from  the  equation  of  state  with  Tm  =  1: 


1 

7(7  -  l)A/r 


(2.23, 


The  S3  and  S4  boundaries  are  treated  in  the  same  manner  as  in  the  first  phase  of  the 

study. 
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Table  2.1.  Phase  2  Boundary  Conditions 


Boundary  Si 

Boundary  S2 

Boundary  S3 

Boundary  S4 

P  =  Po 

P  =  1 

Pz  =  0 

Pr  =  0 

u  =  0 

B23SES1 

u,  =  0 

u  =  0 

v  =  Vr( 2  -  r2)  r  <  1 
t-  =  V'/r  r  >  1 

v  =  V/R 

t>,  =  0 

mm 

tr  =  1 

w  =  1 

w  z  =  0 

t L'r  =  0 

tvz  =  0 

HSBSsSH 

e,  ■=  0 

These  boundary  conditions  were  developed  through  trial  and  error  and  a  more  in¬ 
formative  discussion  on  boundary  conditions  can  be  found  in  Section  5.1.  This  section  is 
intended  only  to  present  the  boundary  conditions  used  in  the  study.  A  summary  of  the 
boundary  conditions  used  in  phase  2  are  in  Table  2.1. 


III.  Solution  Algorithm 


One  of  the  objectives  of  this  study  is  to  be  able  to  determine  if  non  unique  solution.' 
of  the  flow  equations  exist.  To  compute  non-unique  solutions,  a  very  robust  algorithm  is 
needed  to  obtain  both  stable  and  unstable,  steady-state  solutions.  Newton’s  method 
capable  of  computing  non-unique  solutions  for  the  set  of  non-linear,  algebraic  equations 
resulting  from  a  discretization  of  equations  2.11-  2.15. 

Newton’s  method  is  guaranteed  to  converge  if  the  initial  guess  is  within  its  ball  of 
convergence  (19).  The  difficulty  is  finding  an  initial  guess  within  the  ball.  At  higher  values 
of  Reynolds  number  and  vortex  strength,  the  ball  of  convergence  can  contract  making  it 
more  and  more  difficult  to  find  a  suitable  initial  guess.  To  overcome  this  problem  pseudo- 
arclength  continuation  is  used  in  combination  with  Newton's  method.  Pseudo-arclength 
continuation,  like  other  continuation  methods,  uses  information  at  the  current  solution  to 
calculate  a  solution  at  a  parameter  value  close  by.  This  chapter  is  dedicated  to  describing 
these  two  methods. 

3.1  Sewton’s  Method 

Newton’s  method  is  an  iterative  algorithm  that  solves  the  nonlinear  system  of  equa 

tions 

£(x;A)  =  0.  (3.1) 

Given  x\  an  initial  approximation  to  the  solution  vector  x,  and  A  the  free  parameter,  an 
improved  approximation  i,+1  can  be  found  by  solving  the  system  of  equations: 

££U,;A)(x'+1  -s*)  =  -£(*';  A).  (3.2, 


The  solution  to  equation  3.2  is  known  as  a  Newton  iterate.  Successive  Newton  iterates  are 
computed  until  the  norm  of  £  goes  to  zero,  usually  defined  as  machine  precision  or  some 
other  tolerably  small  value.  The  norm  used  is  defined  as: 


II  £11= 


1/2 


ArAzFft  i’+,;A) 


(3.3) 
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F L  is  the  Jacobian  matrix,  which  is  defined  as: 


An  example  calculation  of  some  of  the  elements  of  equation  3. -4  is  shown  in  Appendix  B. 
Newton's  method  is  an  attractive  algorithm  because  of  its  ease  of  programming  and  its 
convergence  rate.  The  method  is  guaranteed  to  converge  quadratically  if  the  Jacobian  ma¬ 
trix  is  non-singular  and  the  initial  guess  is  within  the  ball  of  convergence  (19).  No  singular 
matrices  were  encountered  during  any  of  the  runs  and  pseudo-arclength  continuation  was 
used  to  find  a  suitable  initial  guess. 

3.2  Pscudo-arckngth  Continuation 

It  has  been  proven  that  if  a  solution  £*  is  known  and  £ ^  is  nonsingular,  then  for 
some  range  of  A  about  A*  there  exists  a  unique  solution  path  through  (£*,  A*).  The  proof 
is  outlined  in  (21).  Pseudo-arclength  continuation  (PAC)  is  Keller's  method  of  computing 
solutions  along  the  solution  path  by  using  information  at  (j*,A’)  to  compute  the  next 
solution  point. 

Figure  3.1  is  representative  of  a  solution  path  found  by  plotting  the  norm  of  the 
solution  vector  versus  the  free  parameter  A.  The  PAC  process  is  to  compute  a  tangent 
vector  T  at  the  known  solution  z*.  designated  as  P.  Then,  at  a  distance  d  away,  search 
along  a  line  perpendicular  to  7  for  the  next  solution  point.  This  is  done  by  first  using 
arclength  to  parameterize  the  solution  path  fz  =  z(s),  A  =  A(s).  and  7  =  F(s)  =  0). 
Then  the  tangent  vector  is  computed  by 

j-£U(5);A(s»  =  0.  (3.5) 

ds 

I'sirig  the  chain  rule  equation  3.5  becomes 

£t(>s)j(5)  4-  £x(-s)^(-s)  =  0,  (3.6) 
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Then  equations  3.6  and  3.9  give  the  following  relationships 


AU«) 


±1 

/hTTF 


(3.12) 


and 


i(s)  =  -A(s)<p. 


(3-13) 


The  sign  of  equation  3.12  represents  the  direction  of  the  tangent  vector  and  is  therefore 
indeterminate.  At  the  startup  of  the  continuation  process  this  sign  is  set  depending  on 
which  part  of  the  solution  path  is  to  be  computed. 

From  the  solution  P,  the  tangent  vector  T.  and  the  distance  d ,  the  initial  solution 
vector  Qo  can  be  computed: 


x 

A 


\ 

'  Qo 


(3.1*1  j 


The  solution  Q  lies  on  a  line  perpendicular  to  T  passing  through  Q0 ■  This  condition  can 
be  stated  mathematically  as 


D  =  ip(xg  -  Xj>)  +  (A q  —  Xp)Xp  -  d  (3.15) 

and  can  be  added  to  the  system  of  nonlinear  equations  3.1.  This  new  system  is  then  solved 
by  Newton’s  method  with  Qo  as  the  initial  guess.  Q0  becomes  a  better  and  better  first 
approximation  as  d  gets  smaller.  The  solution  xq  is  obtained  when  F  =  0  and  D-d. 
Another  solution  can  be  computed  by  repeating  the  process. 

For  the  vortex  breakdown  problem,  the  most  informative  graphical  representation  of 
the  solution  path  is  found  by  plotting  the  free  parameter  Re,  V,  or  M  against  the  weighted 
L 2  norm  of  the  change  in  T ( , , ,  from  the  inflow  r0(;): 

(r<i ,j)  -  F0(j))  |  .  (3.16) 

where  T(,  j)  =  r}v,  r  is  referred  to  as  the  “circulation  perturbation”  norm. 
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IV.  Xavier-Stokts  Solutions 


In  Chapter  II.  a  mode!  for  compressible  trailing  vortices  was  developed  utilizing  th>- 
Navier  Stokes  equations.  After  discretization,  the  equations  became  a  set  of  non-linear 
algebraic  equations.  I’ sing  the  method  of  solution  described  in  Chapter  III.  solutions  to 
the  system  of  equations  at  Re  =  200  for  various  vortex  strengths  and  Mach  numbers  were 
computed  and  are  presented  in  this  chapter,  as  well  as  the  solution  paths  computed  by 
continuation  in  \  and  M .  A  grid  sensitivity  study  and  validation  of  the  code  with  previous 
work  are  also  presented. 

4. 1  Gnd  Sensitivity  Study 

In  numerical  studies,  it  becomes  very  important  to  show  the  degree  to  which  the 
finite  difference  solution  is  independent  of  the  grid  used.  Also,  when  using  the  far  field  as  a 
boundary,  the  effect  of  computational  domain  size  on  the  solution  should  also  be  analyzed. 
This  section  consists  of  a  grid  refinement  and  domain  size  study. 

Both  the  studies  were  completed  by  computing  solutions  at  Re  =  200,  V  =  1.0,  and 
M  =  0.15,  since  reversed  flow  occurs  at  this  set  of  parameters.  In  the  grid  refinement 
study  the  maximum  radial  and  axial  distances  were  held  constant  at  R  =  2  and  Z  -  20 
while  the  number  of  nodes  was  varied.  Radial  grid  refinement  was  analyzed  by  holding 
the  nodes  in  the  r  direction  constant  at  /  =  105  and  varying  the  number  of  nodes  in  the 
r  direction  J.  depicted  in  Figure  -4.1.  At  approximately  J  -  27.  the  solution  becomes 
practically  independent  of  increasing  J .  Next.  ;  grid  refinement  was  studied  by  varying 
/  while  holding  J  =  27.  depicted  in  Figure  4.2.  The  solution  points  are  coincidental  for 
/  =  157  and  I  =  209.  The  7  -  105  solution  produces  very  slight  differences  from  the  other 
two  solutions. 

In  the  domain  size  study  the  node  spacing  Ar  =  an(j  =  _iL_  were  J^pt 

approximately  constant  with  R  and  Z  varied.  The  study  of  radial  domain  size  was  difficult 
to  obtain,  because  as  J  is  increased  by  one  the  bandwidth  of  the  system  increases  by  10 
making  the  J  -  40  case  a  computational  limit  for  the  available  resources.  Figure  4.3  shows 
the  comparison  and  although  only  two  solutions  are  shown  (R  =  2  and  R  =  3),  there  is 
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Figure  4.1.  Effect  of  radial-node  spacing  on  centerline  axial  velocity 


2 

Figure  4.2.  Effect  of  axial-node  spacing  on  centerline  axial  velocity 
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Figure  4.3.  Effect  of  radial  domain  length  on  centerline  axial  velocity 


only  a  slight  difference  between  them.  The  R  =  2  solution  has  a  slightly  higher  peak  value 
followed  by  a  slightly  lower  minimum.  The  axial  length  study  was  more  definitive  and  is 
shown  in  Figure  4.4.  At  Z  =  20  and  /  =  105,  the  solution  becomes  independent  of  axial 
domain  length.  In  summary,  the  solution  becomes  practically  independent  of  domain  size 
and  node  spacing  at  R  =  3,  J  -  40,  Z  =  20,  and  I  -  157.  This  would  be  the  grid  to 
choose  if  unlimited  computer  resources  were  available,  but  this  was  not  the  case  for  this 
study.  Instead,  the  standard  grid  used  throughout  the  study  is  R  -  2.  J  =  27,  Z  -  20.  and 
/  =  105.  This  grid  is  a  very  good  compromise  and  has  only  slight  quantitative  differences. 


4-2  Comparison  of  Solutions  for  Low  Mach  Number 

Steady-state  axisymmetric,  trailing  vortices  modeled  with  the  Navier  Stokes  equa¬ 
tions  were  studied  previously  by  Grabowski  (11),  Hafez  et  al.  (12).  Beran  (1),  Hafez  et 
al.  (13).  and  Beran  (2).  Although  each  study  modeled  trailing  vortices  in  essentially  the 
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Figure  4.4.  Effect  of  axial  domain  length  on  centerline  axial  velocity 


same  way,  there  were  significant  differences  in  the  form  of  the  governing  equations,  the 
finite  difference  expressions  used,  and  the  method  of  solution. 

Grabowski  cast  the  governing  equations  in  primitive  variable  form,  similar  to  this 
work,  but  he  used  a  non-uniform  grid  with  refinement  near  the  axis  of  symmetry.  The 
equations  were  integrated  in  time  until  convergence,  utilizing  the  artificial  compressibility 
method  of  Chorin  (4).  Hafez  et  al.  (12)  took  a  different  approach.  He  used  the  stream 
function,  vorticity,  and  circulation  form,  approximating  derivatives  with  upwind  difference> 
at  constant  node  spacing.  An  iterative  relaxation  technique  was  used  to  solve  the  set  of 
equations.  Beran  (1),  Hafez  et  al.  (13),  and  Beran  (2)  also  used  the  stream  function, 
vorticity,  and  circulation  form  with  constant  node  spacing,  but  the  method  of  solving  the 
equations  essentially  follows  Chapter  III  of  this  work. 

All  of  these  studies  were  performed  with  the  assumption  of  M  =  0,  w  hereas  M  =  0  15 
in  the  current  study.  This  Mach  number  was  the  lowest  attainable  due  to  numerical 
difficulties.  This  proved  to  be  a  very  small  portion  of  the  difference  as  seen  through 
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Figure  4.5.  Comparison  of  centerline  axial  velocity  profiles  computed  in  four  investiga 
tions  for  Re  =  200  and  V'  =  1.0 

comparison  with  the  phase  1,  M  =  0  code.  The  qualitative  and  quantitative  behavior  will 
now  be  compared  for  Re  =  200  and  V  =  1.  These  conditions  were  chosen  because  vortex 
breakdown  does  occur  and  data  is  available  from  other  studies  for  these  conditions. 

Figure  4.5  shows  the  comparison  of  these  data  sets.  Solutions  are  in  good  agreement 
until  r  =  2.5.  For  2.5  <  z  <  6  all  five  curves  are  in  what  is  termed  as  phase,  with  variations 
in  peak  value.  Most  notably  Hafez  et  al.  (13)  has  a  much  lower  minimum  value  than  the 
others.  Beran  (2)  and  Grabowski  (11)  are  slightly  lower  in  peak  value.  For  z  >  6  tin- 
present  work,  Grabowski  (11),  and  Hafez  et  al.  (13)  show  excellent  agreement.  Beran  (2) 
returns  to  a  lower  value  of  axial  velocity.  Quantitatively  the  M  =  0.15  and  R  =  3  solution 
and  Grabowski  (11)  are  in  very  good  agreement,  this  is  not  surprising  since  Grabowski  also 
cast  his  equations  in  primitive  variable  form,  whereas  Beran  and  Hafez  used  the  stream 
function,  vorticity,  and  circulation  formulation.  Qualitatively  the  current  study  and  the 
three  previous  studies  all  show  the  same  damped  oscillatory  behavior  with  a  match  in 
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Figure  4.6.  Continuation  in  Mach  number.  Re  =  200  and  V  =  1.0 

phase.  With  this  kind  of  agreement,  the  code  developed  in  this  study  may  be  accepted  as 
a  viable  tool  to  study  the  effects  of  compressibility. 

4-3  Flow  Solutions  at  Various  Mach  Numbers 

In  this  section  solutions  at  varying  Mach  numbers  are  examined.  A  continuation  in 
Mach  number  was  performed  over  the  range  0.2  <  M  <  1.0  with  R(  =  200  and  V'  =  1.0 
In  all  of  these  runs  the  standard  grid  I  =  105,  Z  =  20,  J  =  27,  and  R  =  2  was  used 
Figure  4.6  presents  £p  versus  M  for  the  continuation  run.  As  one  can  see  from  this  figure, 
the  behavior  is  very  regular  as  Mach  number  is  varied.  No  limit  points  were  encountered 
and  the  sign  of  the  determinant  of  the  Jacobian  matrix  was  constant  for  the  complete  run 
suggesting  no  change  in  stability. 

To  analyze  the  effect  of  increasing  Mach  number,  axial  profiles  of  u-  and  p  on  the 
centerline  are  presented  in  Figures  4.7  and  4.8.  In  the  w  profiles,  one  can  see  that  as  Mach 
number  increases  reversed  flow  is  quickly  lost  and  the  breakdown  is  lost  entirely  at  about 
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Figure  4.9.  v  profiles  on  the  centerline  for  various  Mach  numbers.  Re  =  200  and  V  =  1  0 


M  =  0.5.  Thus,  by  increasing  Mach  number,  breakdown  is  effectively  delayed  for  V  held 
constant. 

To  explain  this  trend  a  look  at  the  relevant  terms  in  the  radial  momentum  equation 
proves  fruitful.  Beran  (2)  found  that  with  the  assumption  of  a  quasi -cylindrical  flow,  the 
relevant  terms  in  the  incompressible  flow  radial  momentum  equation  are 


t’2  dp 

7  *  Tr' 


(4.1  i 


Similarly  for  compressible  flow,  the  radial  momentum  equation  becomes 

pv2  dp 
—  * 
r  dr 


(4.2.i 


The  assumption  of  quasi-cylindrical  flow  is  not  a  perfect  assumption  but  certainly  retains 
the  largest  terms  of  the  radial  momentum  equation  and  becomes  more  applicable  as  M 
gets  larger,  as  can  be  observed  in  Figure  4.7.  Integrating  equation  4.2  with  respect  to  r 
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across  the  domain  obtains  the  relationship 


but  is  equal  to  zero  if  R  is  large  enough.  So,  equation  4.3  becomes 


Po=  ~ 


(4.3) 


(4.4) 


To  obtain  the  axial  pressure  gradient  on  the  centerline,  take  the  partial  derivative  of 
equation  4.4  with  respect  to  z,  obtaining 


dp0 

dz 


(4.5.1 


Beran  (2)  found  that  for  z  between  the  inflow  and  the  axial  position  of  minimum  axial 
velocity,  an  adverse  pressure  gradient  caused  a  stagnation  of  the  flow.  The  terms  J,  r2, 
and  2pv  are  positive  definite,  so  the  terms  ||  and  become  the  terms  that  can  cause 
a  favorable  pressure  gradient  or  an  adverse  pressure  gradient  for  a  compressible  flow. 
Figure  4.9  shows  the  axial  profiles  of  v  at  R  =  1  and  for  M  =  0.2  and  0.4.  Both  the 
M  =  0.2  and  0.4  profiles  have  a  negative  |j,  with  its  magnitude  getting  smaller  as  M 
is  increased.  The  second  part  of  equation  4.5  then  provides  an  adverse  pressure  gradient 
whose  magnitude  gets  smaller  as  M  is  increased.  The  axial  density  profiles  in  figure  4.8 
provide  information  on  the  first  term  in  equation  4.5.  The  term  ||  is  positive  and  gets 
larger  as  M  is  increased.  This  again  decays  the  adverse  pressure  gradient  by  providing  a 
negative  |*  contribution.  In  summary,  as  A/  is  increased  the  adverse  pressure  gradient,  by 
which  the  axial  velocity  on  the  centerline  is  stagnated,  is  reduced. 

It  is  also  interesting  to  note  that  as  A/  is  increased  at  regular  intervals  the  density 
profiles  change  in  like  manner.  This  relationship  explains  the  near  linear  behavior  of 
Figure  4.6.  A  solution  was  obtained  at  Af  =  1.0  without  difficulty,  even  though  the 
non  conservative  form  of  the  equations  was  used.  No  spikes  in  the  velocity  profiles  were 
observed.  The  flow  solution  is  composed  of  weak  shock  waves.  The  assumption  of  a  linear 
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Figure  4.10.  Continuation  in  V  for  Re  =  200  and  M  =  0.2 

viscosity  temperature  relationship  produces  a  maximum  error  of  5.0%  at  this  Mach  number 
however. 

4-4  V  Continuation  at  Various  Mach  Numbers 

One  of  the  important  discoveries  of  Beran  (2)  was  the  existence  of  non-unique  solu¬ 
tions  for  a  range  of  vortex  strengths  and  Re  =  200.  One  of  the  objectives  of  this  work  is  to 
see  what  happens  to  the  non-unique  solutions  as  Mach  number  is  increased.  An  important 
first  step  is  to  reproduce  the  non-unique  behavior  as  M  —  0. 

Figure  4.10  shows  a  continuation  run  at  M  =  0.2,  Re  =  200,  and  0.5  <  V'  <  3.0.  The 
non-unique  behavior  does  indeed  exist  at  this  low  Mach  number.  A  change  in  sign  of  the 
Jacobian  matrix  determinant  was  found  at  the  limit  point,  which  generally  is  associated 
with  a  change  in  stability.  At  V  =  0.9,  the  continuation  process  broke  down.  This  was 
probably  due  to  a  second  limit  point  with  too  large  a  value  of  d  to  resolve.  In  an  attempt 
to  fill  in  this  curve  the  M  =  0.4  and  V  =  1.5  solution  was  used  to  obtain  an  M  =  0.2 
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V 

Figure  4.13.  Continuation  in  V  for  Re  =  200  and  M  =  0.2, 0.4.  and  0.6 


and  V  =  1.5  solution  through  continuation  in  M .  This  solution  was  then  used  to  begin 
a  continuation  run  down  from  V  =  1.5.  Another  limit  point  was  found  near  the  previous 
one  at  V  =  1.0589.  The  sign  of  the  determinant  changed  through  the  second  limit  point. 
The  curves  were  not  able  to  be  connected  because  of  the  breakdown  of  the  continuation 
process.  Continuation  runs  were  then  made  for  M  =  0.4  and  M  =  0.6.  These  curves 
are  shown  in  Figures  4.11  and  4.12.  Non  unique  solutions  were  not  evident  in  these  runs 
and  no  change  in  sign  of  the  determinant  was  found.  Both  solution  paths  were  taken  to 
extremely  high  vortex  strengths  to  ensure  the  absence  of  non-unique  solutions. 

Figure  4.13  is  a  composite  plot  of  the  three  continuation  runs  in  the  region  of  non 
unique  solutions.  The  three  solutions  are  practically  coincidental  up  to  a  vortex  strength 
of  approximately  0.75.  At  V  =  0.75  the  three  curves  separate.  The  M  =  0.2  path  makes 
an  abrupt  slope  change  and  then  runs  a  parallel  course  with  the  other  two  solution  paths 
until  V'  =  1.0587.  At  V  =  1.0587  a  limit  point  is  encountered.  The  .Vf  =  0.4  solution 
exhibits  some  of  the  slope  change  at  V'  =  0.75  as  well,  but  not  as  abrupt  a  change  for 
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Figure  4.14.  w  profiles  on  the  centerline  for  various  V’s  at  Re  =  200  and  M  =  0.8 


M  =  0.2.  The  M  -  0.6  curve  makes  an  even  smaller  slope  change  and  then  runs  a  parallel 
course  to  the  M  =  0.4  path. 

Figures  4.14  and  4.15  show  the  centerline  axial  velocity  and  density  profiles  of  four 
increasing  vortex  strengths  at  M  =  0.8.  These  solutions  are  presented  to  show  the  effects 
of  increasing  V'  at  a  high  Mach  number.  As  V  was  increased  the  minimum  value  of  u 
moves  towards  the  inflow  boundary  until  a  certain  value  of  Y  is  reached.  For  M  =  0.8 
this  value  is  between  1.0  and  1.2.  At  this  point  increasing  V7  compresses  the  fluid  on  the 
centerline  with  slight  increases  in  the  minimum  u\  This  increase  in  minimum  u-  was  due 
to  an  enforcement  of  conservation  of  mass  as  density  becomes  lower  at  this  axial  position. 
There  was  no  reverse  flow  observed  for  this  Mach  number. 
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V'.  Conclusions  and  Recommendations 


This  chapter  discusses  the  results  presented  in  Chapter  4  and  the  difficulties  encoun¬ 
tered  in  finding  the  proper  boundary  conditions.  A  summary  of  the  effects  of  compressibil- 
itv  on  the  vortex  breakdown  and  a  presentation  of  recommendations  for  follow  on  study 
are  given. 

5. 1  Solution  Sensitivity  to  Boundary  Conditions 

Determining  the  boundary  conditions  for  the  compressible  vortex  was  not  a  trivial 
task.  Unlike  mam  fluid  dynamics  problems,  the  inflow  conditions  are  at  a  station  where 
flow  conditions  are  not  necessarily  known.  The  conditions  on  u.  v,  and  w  were  taken  from 
Beran's  incompressible  study  (2).  Conditions  that  specify  the  inflow  density  and  inter¬ 
nal  energy  profiles  become  very  difficult  to  obtain.  The  first  set  of  conditions  tried  were 
columnar  flow  conditions,  pz  =  =  et  -  0.  Since  the  typical  solution  to  the  incompress¬ 

ible  flow  equations  is  damped  oscillatory  motion,  specification  of  a  locally  columnar  flow 
should  have  set  the  conditions  at  the  peak  of  one  of  the  flow  oscillations.  The  result  of 
these  conditions  was  immediate  divergence  even  at  very  low  values  of  V .  This  suggests  an 
overconstraint  in  the  conditions. 

To  alleviate  this  problem,  the  conditions  on  p  and  e  were  changed.  A  density  profile 
was  obtained  by  solving  the  columnar  flow  equations  and  used  as  a  Dirichlet  condition 
on  p.  Then  the  u  momentum  equation  was  evaluated  on  the  inflow  boundary  with  ax¬ 
ial  derivatives  approximated  by  first-order  difference  operators.  This  was  necessary  since 
second-order  operators  would  increase  the  bandwidth  by  5.  which  was  unacceptable  compu¬ 
tationally.  These  conditions  worked  for  lower  values  of  V,  but  at  V'  =  0.8,  large-amplitude 
numerical  noise  appeared  at  the  inflow  boundary  and  grew  worse  as  V  was  increased. 
Satisfying  the  u  momentum  equation  at  the  inflow  boundary  as  the  condition  on  p,  and 
satisfying  the  r  momentum  equation  as  the  condition  on  e,  was  also  tried  with  similar 
results. 

It  seemed  as  though  the  lack  of  information  on  c  was  causing  the  difficulty,  so  the 
next  set  of  conditions  tried  were  conditions  on  p,  u,  v,  and  w  alone,  allowing  e  to  be 
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set  by  the  governing  equations.  Specifically,  these  conditions  were  p  =  p o  obtained  from 
the  columnar  flow  equations,  u  =  0,  v  set  by  Figure  2.2,  u-  =  1.  and  tu,  =  0.  These 
conditions  worked  well  for  Re  =  200  and  the  range  of  V’"s  physically  attainable  and  were 
the  conditions  used  throughout  this  work.  The  iterative  solution  process  diverges  for  higher 
Reynolds  numbers  due  to  the  boundary  conditions  though,  hence  the  reason  Re  =  200  was 
used  throughout  the  study.  A  more  robust  set  of  conditions  need  to  be  determined.  Either 
experimental  p  and  e  profiles,  or  an  extension  of  the  quasi-cylindrical  equations  found  in 
Beran  (2)  to  compressible  flow  should  be  utilized.  In  the  absence  of  vortex  breakdown,  the 
quasi-cylindrical  equations  approximate  trailing  vortices  at  high  Reynolds  numbers  very 
well.  For  this  reason,  they  should  provide  a  good  set  of  initial  conditions  for  p  or  e. 

5.2  Compressibility  Trends 

There  were  some  very  interesting  trends  observed  in  this  compressibility  study.  Com¬ 
pressibility  has  a  favorable  effect,  from  the  viewpoint  of  an  aerodynamicist,  on  vortex  break¬ 
down.  This  can  be  observed  in  Section  4.3;  as  Mach  number  was  increased,  breakdown 
was  quickly  lost.  This  was  probably  due  to  compression  of  the  fluid  reducing  the  adverse 
pressure  gradient  on  the  centerline,  destroying  the  the  vortex  breakdown  phenomenon. 

In  Section  4.4,  the  effect  of  varying  vortex  strength  for  a  constant  Re  and  M  was 
analyzed.  When  flow  solutions  at  M  =  0.8  were  computed,  reversed  flow  was  not  observed, 
even  for  extremely  large  values  of  V. 

Beran 's  incompressible  study  (2)  showed  the  occurrence  of  non-unique  solutions  for 
increasing  V  at  Re  =  200,  and  at  realistic  vortex  strengths.  As  Mach  number  was  increased, 
these  non-unique  solutions  disappeared.  This  is  desirable  because  multiple  solutions  cause 
problems  with  time-integration  schemes,  since  the  observed  steady  state  solution  becomes 
dependent  on  the  initial  flow  state. 

In  summary,  as  Mach  number  was  increased  for  Re  =  200  and  V  =  1.0,  breakdown 
was  lost  at  about  M  =  0.5.  As  vortex  strength  was  increased  for  A/  =  0.8,  the  minimum 
value  of  w  moved  towards  the  inflow  boundary,  as  was  reported  by  Beran  (2).  In  Beran's 
incompressible  study  (2)  as  V'  was  increased  other  local  minima  were  observed.  In  the 
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compressible  study,  when  V  grew  large  enough  the  minimum  w  stopped  moving  toward 
the  inflow  ,  but  instead  of  producing  more  local  minima  as  in  incompressible  flow,  the  fluid 
becomes  more  compressed. 

5.3  Recommendations  for  Follow-on  Study 

During  the  course  of  this  study  original  goals  of  mapping  out  the  solution  space  for 
various  values  of  Re,  V ,  and  M  had  to  be  reduced  to  accomplishing  a  compressibility  study 
at  low  Reynolds  number  for  various  vortex  strengths.  This  was  due  to  time  constraints  and 
the  difficulties  encountered  with  the  boundary  conditions.  One  can  see  that  a  logical  follow  - 
on  study  would  be  to  derive  a  set  of  robust  boundary  conditions  that  can  accommodate 
larger  Reynolds  numbers,  as  well  as  various  Mach  numbers  and  vortex  strengths.  This 
follow-on  study  would  benefit  time  integration  studies  as  well,  since  time-integration  codes 
may  converge  on  a  solution  even  if  the  wrong  set  of  boundary  conditions  are  used.  The 
compressible  quasi-cylindrical  equations  mentioned  previously  are  a  promising  starting 
point. 

Once  these  boundary  conditions  are  determined,  simulations  at  higher  Reynolds  num¬ 
bers  should  be  run  to  determine  whether  non-unique  solutions  can  be  found.  Also,  flow 
structures  should  be  studied  to  determine  if  the  same  trends  are  found  for  high  Reynolds 
numbers  as  was  found  for  Re  =  200. 

Finally,  a  compressible  time-integration  code  should  be  developed  to  determine  the 
stability  of  the  different  branches  in  the  continuation  curves.  With  the  current  study  and 
the  follow  on  topics  providing  the  compressible  flow  information  on  vortex  breakdown,  one 
would  hope  that  future  investigators  would  be  able  to  put  the  experimental,  analytical, 
and  numerical  studies  together  to  form  a  better  understanding  of  the  vortex  breakdown 
problem. 
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Appendix  A.  Derivation  of  the  Governing  Equations 


The  derivation  of  the  governing  equations  was  performed  by  utilizing  the  nondimen- 
sionai  genera]  curvilinear  coordinate  Navier  Stokes  equations  found  in  reference  (34).  The 
general  continuity  equation  is: 


dp  1  d 
dt  +  hih2h3  dx} 


[^-hxhihzpiij 


) 


=  0 


( A .  1 ) 


where  the  j  subscript  denotes  summation.  Cylindrical  coordinates  were  used  in  the  model 
In  cylindrical  coordinates  x^  =  r.  x2  =  9,  and  x3  -  z  with  velocity  components  u;  equal 
to  u.  v,  and  tv.  The  metrics  h}  are  as  follows:  hj  =  1,  h2  =  r,  and  h3  =  1.  Equation  A.l 
becomes: 


dp  1 

di^r 


d 


d 


^(rp*)+-{pv)+-(rpuO 


=  0 


(A. 2) 


Two  of  the  main  assumptions  of  this  thesis  are  ^  =  0  (steady  state)  and  ^  =  0  (ax- 
isy  mmetric).  With  these  assumptions  applied  to  Equation  A. 2  the  cylindrical  axisymmetric 
continuity  equation  becomes: 


d  pit  d 

_,p„)  + _  +  ,  =  0 


(A.3) 


The  general  curvilinear  coordinate  momentum  equation  in  the  xj  direction  can  be 
written: 


dt 


(pu i)  + 


dh2 


(A. 4) 


with  the  other  two  momentum  equations  following  by  cyclic  permutation  and  where  Tx: 
are  the  components  of  the  tensor  T  : 


TXj  —  pu,u:  -f  pi tj  — 


(A. 5) 
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(defV)u  =  ur 

(A. 11) 

(defV)n  =  - 

r 

(A. 12) 

(de/V)  33  =  w. 

(A. 13) 

.......  1  /  #U’\  1 

(^/V)31  =  -  +  =  _(u.,  +uv) 

(A  14. 

r31  =  2*2  (  U2  +  UV  )  =  *i  ( u,  +  uv  ) 

(A. 15) 

=  “^  (ur  +  “  +  U',  j  +  2*iU, 

(A. 16) 

2  /  u  \  u 

^22  =  ~3^  (  Ur  +  7  +  U’1  )  +  2^“ 

(A. 17; 
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After  making  the  assumption  of  steady  state  and  axisymmetry  Equation  A. 4  becomes: 


-|-<r^i)+  -|-(f-/3I)-  —  =  0 
r  dr  r  dz  r 


using  the  relationships  A. 10-  A. 17  the  relevant  T  terms  become: 


[  A.  l*i 


Tn  -  puu  +  p  -  — 


(^Ur  +  ^  +  2/itir 


/■3,  =  puw  -  ~[p(ut  +  Wr  )] 


'  2  /  u  \  A  u 

-~H  ^uT  -f  -  +  tr.  j  t  2p- 

Substituting  these  terms  into  A.  18  the  u  momentum  equation  becomes: 


Tu  -  P ft’  +  P  -  7^ 


(A. 19) 
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Subtracting  continuity  and  combining  terms: 

pv2 

puuT - —  +  pwu,  +  pT  ~ 

-^-[pr(4/3ur  -  2/3—  -  2/3u>,)  +  pt(u,  +  uv) 
He  t 

+p(4/3urr  +  4/3 -  4/3  — j  +  +  l/3uv*)] 

In  Section  2.1  the  equation  of  state  was  given  to  relate  p,  p.  and  e: 


(A. 22) 


p  =  h  -  l)pc  (A. 23) 

The  linearized  viscosity-temperature  relationship  was  also  given  to  relate  p  and  e: 

p  =  cjf  +  c2  (A. 24  ) 

where  cj  and  c2  are  defined  by  Equations  2.16  and  2.17.  Substituting  these  relationships 
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into  A. 23  the  u  momentum  equation  is  obtained: 


Pv 

PUUr - —  +  P U-U;  +  ("1  -  1  )(pe) r  - 

4-[cif,(4/3u,  -  2/3—  -  2/3u', )  -*■  c,e,(ut  +  uv  )  (  A  2o  , 

fit  T 

+  (cif  4-  Cj)(4/3urr  4-  4/3 —  —  4/3  — j  4-  u,,  4-  1/3 txv * )] 


The  v  momentum  equation  found  by  cyclic  permutation  and  the  application  of  pre¬ 
viously  mentioned  assumptions  is: 


ia  r  1  d  ^ 

-x-(rj^]2)+  -  r-^32 )  +  - 

r  or  r  dz  r 


=  0 


(A. 26 ) 


For  the  r  momentum  equation  more  of  the  ( defV)l}  and  rtJ  terms  must  be  defined: 


(defV) „  =  (de/V)2l  =  l-  [vr  -  ^ 

(A. 27) 

(de/V)32  =  ~(vt) 

(A. 28) 

(  v\ 

rl  2  =  r21  =  P  1  Vr  -  -  1 

(A. 29) 

r32  =  P(V,) 

(A. 30) 

Using  A. 27-  A. 30  the  relevant  Tx]  terms  become: 

^12  =  ■/"< 21  =  P™  -  ^  ^  ^tv  - 

(A  31 ) 

J32  =  PfU’  ~  Ye 

(A. 32) 

Substituting  these  terms  into  A. 26  the  equation  becomes: 
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Subtracting  continuity  and  combining  terms: 


put' 

pu  tv  +  - +  pU'V,  = 

r 

1  v 

[Pr(  Vr  ~  —  )  "t"  Pi  (  *■'»  ) 

Re  r 

,  VT  f 

+p(v„  + - 5  +  I’tl 

r  rl 


( A. 33 ) 


Using  Equation  A. 24  the  r  momentum  equation  becomes: 


puv 

pu  tv  +  - +  pwi  - 

T 

1  t’ 

—  [cjer(iv  -  -)  +  c,ef(Vj) 

/if  r 

+  (Clf  +  C2)(l'rr  +  — - j  +  t’,j)] 

t  r * 


( A. 34  j 


The  u>  momentum  equation  found  by  cyclic  permutation  and  the  application  of  pre¬ 
viously  mentioned  assumptions  is: 


1  d  ,  _  4  1  6  .  _  .  „ 

757(r-F,3> +  = 0 


(A. 35) 


For  the  u-  momentum  equation  more  of  the  (defV)tJ  and  rtJ  terms  must  be  defined: 


(defV)  13  =  ~(wr  +  uf) 


(A. 36) 


^3  =  M(u-r  +  u,) 

/  2  2  u  4  \ 

r33  =  ^-5«,-57+3«-,j 

Using  A. 36-  A. 38  the  relevant  TX]  terms  become: 


(A. 37) 
(A.  38) 


13  =  pw  -  +  «*)] 


«<.»»  + P -  i  [*(-!«--  57  +  5^)] 


(A. 39) 

(A. 40) 
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Substituting  these  terms  into  A. 35  the  equation  becomes: 


d 

1  , 

1  1  1 

dr 

putr  -  —  (puy  +  pu,) 
Re  J 

+  -  ptltt’  —  —  (pU'r  4-  fiUr  ) 

r  Re 

1 

(  2 

2  u  4  \1 

r  -  Te 

-  Vt  + 

Subtracting  continuity  and  combining  terms: 


putty  +  pww,  +  px  = 

-^-[pr(u’T  +  u,)  +  p,(-2/3ur  -  2/3—  +  4/3u’z)  (A. 41) 

He  t 

+m(1/3u,j  +  l/3y  +  u-Vr  +  y  +  4/3u'jj  )] 


Using  Equation  A. 24  the  tr  momentum  equation  becomes: 


putty  +  pww,  +  (7  -  l)(pc)x  = 

[c  j  ( tty  +  u,)  +  Cie,(-2/3ur  -  2/3-  +  4/3u-x )  (A. 42) 

Re  r 

+  (Cje  +  Cj)(  l/3ur,  +  1/3^  +  tty,  +  +  4/3u>zl)J 


The  non-dimensional  genera]  curvilinear  coordinate  energy  equation  can  be  written: 


dtipE)+  h^^dz 


;R 


,  r  .  \  1  -jk  1  de 


=  0  (A  43) 


where  £  =  e  +  i(u2  +  v2  +  tr2).  Now  applying  the  above  assumptions  and  metrics  A. 43 
becomes: 


d_ 

dr 


1  ■>  2  2  1  7  k  de 

pue  +  -(pu  +  pur  +  puw* )  +  pu  -  —  (rnur,2r  +  r13ii')  - 

d  f  1  2  7  3a  1  .  -)k  de 

+  d~-  \ptu' +  2(pu  w  +  pv  w  +  pU'  ^  +  pu'  ~  ^T3lU  +  T3iV  +  r*»u') "  rTPt O' 


(A. 44 
•jk  de 


=  0 


The  stress  tensor  terms  in  Equation  A. 45  have  ail  been  defined  above.  Applying  these 
terms  and  using  A. 23  A. 45  becomes: 
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■>!  —  +  (pue)r  4-  (pri'e),)  +  y[ti3  +  tit-2  +  uu>2]  4  y|u2u>4-  t,2*-  +  u-3] 

+  ^I(^3)r  +  («t’2)r  +  (UWJ)T  +  (ti2ttf),  +  (v2U'),  +  (ti>3),  +  —  +  —  +  ^-1  = 

1  rrr 

1  ,  U2 

—  [cjer(4/3uur  -  2/3—  -  2/3uu-,  4-  t-t-r  -  —  +  u*tr  4-  truv) 

+  C\eI(uui  +  uuv  4-  t-t-,  -  2/3t ltw  -  2/3—  4-  4/3u-u,)  (A. 45) 

+  (cie  4-  c2)(4/3u2  +  4/3uurr  -  4/3tirti>,  4-  1/3 utc„  4-  l/3uriu'  4-  t-2 

VVT 

4-  vvTT - —4-  2tijtur  4-  u>,  4-  tuu)rr  4-  u2  4-  uu,z  4-  r2  4-  vvtz 

+  1/3^  -  4/3^  +  4/3u-2  +  4/3u«u-„  -  2/3y  -  ~  + 

+  (cic  +  C2)(crr  4-  e!2  4-  y )  4-  Cje2  4-  c2  ] 
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Appendix  B.  Derivation  of  Representative  Analytic  Jacobian  Elements 

The  purpose  of  this  appendix  is  to  set  up  Newton’s  method  for  a  representative 
equation.  The  process  of  going  from  a  differentia]  equation  to  an  algebraic  equation  and 
then  computing  the  associated  analytic  Jacobian  elements  is  described  in  detail. 

Due  to  its  simplicity  the  continuity  equation  will  be  used  to  demonstrate  the  process. 
The  differential  form  of  the  continuity  equation  is: 

pu 

(Pu)r  +  —  +  (PU>)X  =  0  (B.l) 

The  first  step  is  to  put  the  equation  in  non-conservative  form: 

pu 

pru  +  pur  +  —  +  ptw  +  pii't  =  0  (B.2) 

In  the  discretization  process  all  partial  derivatives  in  the  interior  equations  are  replaced 
with  discrete  operators.  In  a  2-D  implicit  scheme  with  J  nodes  in  the  r  direction,  I  nodes 
in  the  z  direction,  and  J  <  I  a  column  by  column  numbering  system  produces  the  smallest 
bandwidth  of  the  banded  system  to  be  solved.  The  stencil  used  is  found  in  Figure  B.l. 
With  this  stencil,  the  second  order  central  difference  operator  is: 


£  «kt  -  «fc3 

6,u  =  - 

2  Az 


Applying  the  6r  and  Sz  operators  to  each  of  the  terms  in  Equation  B.2  it  becomes: 


Pk2  ~  PkA  ,  2  ~  ,  Pk^k  ,  Pk  1  “  Pk3  ,  U  fcl  ~  „  ,  , 

nr|  “* +l,k  +  — +  +l,t  Hit-  = 0  lB  4! 


This  is  the  non-linear  algebraic  continuity  equation  to  be  solved  by  Newton’s  method  and 
will  be  referred  to  as  F. 

The  next  step  is  to  compute  the  analytic  Jacobian  elements  of  B.4.  This  is  done  by 
taking  partial  derivatives  of  B.4  with  respect  to  each  variable  at  each  node.  The  Jacobian 
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(B.ll) 


dF  _  1 

dun  2A  rPk 

dF  _  1 

duk  4  2A  rPk 

and  the  elements  associated  with  tr  are: 


(  B .  1 2  i 


dF _ 1_ 

dwki  2A  z^k 


( B .  1 3 ) 


dF 


1 


■Pk 


[B.14) 


dwk  3  2Ar 

This  process  is  repeated  for  the  momentum  equations  and  the  energy  equation,  and  then 
Newton’s  method  is  used  to  solve: 


£*Ai  =  -£ 


(B.15 ) 


where  £ ^  is  the  Jacobian  matrix. 
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